closs_tbsvm <- function(X, y, lambda1 = 1, lambda2 = lambda1, 
                        gamma1 = 1e-7, gamma2 = gamma1, 
                        kernel = c("linear", "rbf", "poly"),
                        gamma = 1 / ncol(X), degree = 3, coef0 = 0,
                        sigma1 = 1, sigma2 = sigma1,
                        eps.hq = 1e-2, hq.steps = 10,
                        fit_intercept = TRUE,
                        randx = 1) {
  closs_tbsvm_dual_solver <- function(KernelX, idx, lambda1, lambda2, 
                                      gamma1, gamma2,
                                      sigma1, sigma2,
                                      eps.hq, hq.steps) {
    update_omega <- function(f, sigma_sq) {
      return(as.vector(exp(-f^2/(2*sigma_sq))))
    }
    
    lambda1_ <- lambda1/(2*sigma1^2)
    lambda2_ <- lambda2/(2*sigma2^2)
    
    X1 <- KernelX[idx, ]
    X2 <- KernelX[-idx, ]
    
    xn <- nrow(KernelX)
    xp <- ncol(KernelX)
    X1n <- nrow(X1)
    X2n <- xn - X1n
    
    Omegavec1 <- rep(1, X2n)
    Omegavec2 <- rep(1, X1n)
    
    gamma1I <- diag(gamma1, xp)
    gamma2I <- diag(gamma2, xp)
    
    e1 <- matrix(1, X1n)
    e2 <- matrix(1, X2n)
    
    u01 <- matrix(0, xp)
    u02 <- matrix(0, xp)
    
    GramX1 <- t(X1) %*% X1
    GramX2 <- t(X2) %*% X2
    for (i in 1:hq.steps) {
      X2omegalam_T <- t(X2*(Omegavec1*lambda1_))
      u11 <- -cholsolve((GramX1 + X2omegalam_T %*% X2 + gamma1I), 
                        X2omegalam_T %*% e2)
      f1 <- X2 %*% u11
      Omegavec1 <- update_omega(1 + f1, sigma1^2)
      if (norm(u11 - u01, type = "2") < eps.hq) {
        break
      } else {
        u01 <- u11
      }
    }
    for (i in 1:hq.steps) {
      X1omegalam_T <- t(X1*(Omegavec2*lambda2_))
      u12 <- cholsolve((GramX2 + X1omegalam_T %*% X1 + gamma2I), 
                       X1omegalam_T %*% e1)
      f2 <- X1 %*% u12
      Omegavec2 <- update_omega(1 - f2, sigma2^2)
      if (norm(u12 - u02, type = "2") < eps.hq) {
        break
      } else {
        u02 <- u12
      }
    }
    
    BaseDualClossTBSVMClassifier <- list("coef1" = as.matrix(u11),
                                         "coef2" = as.matrix(u12))
    return(BaseDualClossTBSVMClassifier)
  }
  X <- as.matrix(X)
  y <- as.matrix(y)
  class_set <- sort(unique(y))
  idx <- which(y == class_set[1])
  y[idx] <- -1
  y[-idx] <- 1
  y <- as.matrix(as.numeric(y))
  if (length(class_set) > 2) {
    stop("The number of class should less 2!")
  }
  kernel <- match.arg(kernel)
  solver <- "dual"
  kso <- SupportVectorLab:::kernel_select_option(X, kernel, "primal", randx,
                                                 gamma, degree, coef0)
  KernelX <- kso$KernelX
  Kw <- kso$KernelX[kso$sample_idx, ]
  X <- kso$X
  if (fit_intercept == TRUE) {
    KernelX <- cbind(KernelX, 1)
  }
  if (solver == "dual") {
    solver.res <- closs_tbsvm_dual_solver(KernelX, idx, lambda1, lambda2, 
                                          gamma1, gamma2,
                                          sigma1, sigma2,
                                          eps.hq, hq.steps)
  }
  TSVMClassifier <- list("X" = X, "y" = y, "class_set" = class_set,
                         "lambda1" = lambda1, "lambda2" = lambda2, "kernel" = kernel,
                         "gamma" = gamma, "degree" = degree, "coef0" = coef0,
                         "solver" = solver, "coef1" = solver.res$coef1,
                         "coef2" = solver.res$coef2,
                         "fit_intercept" = fit_intercept,
                         "Kw" = Kw)
  class(TSVMClassifier) <- "TSVMClassifier"
  return(TSVMClassifier)
}